GARCH(1,1) Models (Volatility Models)#
Goals#
Precisely define conditional heteroskedasticity and connect it to volatility clustering.
Explain why volatility (conditional variance) is modeled rather than the mean for many return series.
Derive and break down the GARCH(1,1) variance recursion in LaTeX.
Implement simulation, conditional variance recursion, and variance forecasts in low-level NumPy.
Visualize volatility clustering and variance forecasts using Plotly.
Conditional heteroskedasticity (quick recap)#
A return innovation \(\varepsilon_t\) is conditionally heteroskedastic if its conditional variance depends on past information:
In finance, \(h_t\) is interpreted as volatility squared (conditional variance).
Why model variance (volatility) instead of the mean?#
Mean returns are often weakly predictable at short horizons; volatility is typically persistent.
Many tasks require a forecast of risk (uncertainty), not just expected return.
Common financial use cases:
VaR / Expected Shortfall and stress testing
Volatility forecasting for portfolio risk budgets and leverage control
Options / derivatives context (realized vs implied volatility)
Position sizing (e.g., target volatility strategies)
GARCH(1,1) model (LaTeX)#
GARCH extends ARCH by adding lagged conditional variance (persistence).
Variance equation breakdown#
\(\alpha\): how strongly “new information” (yesterday’s squared shock) moves volatility.
\(\beta\): how persistent volatility is (how slowly it decays).
Common parameter constraints#
\(\omega > 0\), \(\alpha\ge 0\), \(\beta\ge 0\)
Second-order stationarity (finite unconditional variance): \(\alpha + \beta < 1\)
Under \(\alpha + \beta < 1\), the long-run variance is: $\( \bar h = \mathbb{E}[h_t] = \frac{\omega}{1-\alpha-\beta}. \)$
Volatility clustering + market-shock intuition#
A big market shock at time \(t\) means \(|\varepsilon_t|\) is large.
That makes \(\varepsilon_t^2\) large, which pushes up \(h_{t+1}\).
If \(\beta\) is large, elevated variance persists for many periods, producing clusters of large moves.
Variance forecasts (multi-step)#
Given \(\mathcal{F}_t\), the one-step-ahead forecast is directly:
For \(k\ge 2\), we use \(\mathbb{E}[\varepsilon_{t+k-1}^2\mid\mathcal{F}_t] = h_{t+k-1\mid t}\), giving the recursion: $\( h_{t+k\mid t} = \omega + (\alpha+\beta)\,h_{t+k-1\mid t}. \)$
Closed form (mean reversion to \(\bar h\)): $\( h_{t+k\mid t} = \bar h + (\alpha+\beta)^{k-1}\,(h_{t+1\mid t} - \bar h). \)$
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import os
import plotly.io as pio
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
pio.templates.default = "plotly_white"
np.set_printoptions(precision=4, suppress=True)
def simulate_garch11(
T: int,
omega: float,
alpha: float,
beta: float,
mu: float = 0.0,
burn: int = 500,
seed: int = 7,
):
"""Simulate a GARCH(1,1) process with Gaussian innovations."""
if omega <= 0:
raise ValueError("omega must be > 0")
if alpha < 0 or beta < 0:
raise ValueError("alpha and beta must be >= 0")
if alpha + beta >= 1:
raise ValueError("Need alpha + beta < 1 for finite unconditional variance in this demo")
h_bar = omega / (1.0 - alpha - beta)
n = int(T + burn)
rng = np.random.default_rng(seed)
z = rng.standard_normal(n)
eps = np.zeros(n)
h = np.full(n, h_bar)
eps[0] = np.sqrt(h_bar) * z[0]
for t in range(1, n):
h[t] = omega + alpha * (eps[t - 1] ** 2) + beta * h[t - 1]
eps[t] = np.sqrt(h[t]) * z[t]
eps = eps[burn:]
h = h[burn:]
r = mu + eps
return r, eps, h
def garch11_conditional_variance(eps, omega: float, alpha: float, beta: float, initial_variance: float | None = None):
eps = np.asarray(eps, dtype=float)
if initial_variance is None:
if alpha + beta >= 1:
raise ValueError("Need alpha + beta < 1 to use the unconditional variance as initialization")
initial_variance = omega / (1.0 - alpha - beta)
h = np.full(eps.size, float(initial_variance))
for t in range(1, eps.size):
h[t] = omega + alpha * (eps[t - 1] ** 2) + beta * h[t - 1]
return h
def garch11_forecast_variance(eps_last: float, h_last: float, omega: float, alpha: float, beta: float, horizon: int):
"""Forecast h_{t+k|t} for k=1..horizon."""
horizon = int(horizon)
if horizon < 1:
return np.array([], dtype=float)
h_fore = np.zeros(horizon, dtype=float)
h_fore[0] = omega + alpha * (eps_last**2) + beta * h_last
for k in range(1, horizon):
h_fore[k] = omega + (alpha + beta) * h_fore[k - 1]
return h_fore
# Example: simulate a persistent GARCH(1,1) to show volatility clustering
T = 2500
mu = 0.0
omega = 0.02
alpha = 0.08
beta = 0.90 # alpha+beta close to 1 => strong persistence
r, eps, h = simulate_garch11(T=T, omega=omega, alpha=alpha, beta=beta, mu=mu, seed=123)
view = 800
df = pd.DataFrame(
{
"t": np.arange(T),
"return": r,
"eps_sq": eps**2,
"h": h,
"sigma": np.sqrt(h),
}
).tail(view)
fig = make_subplots(
rows=2,
cols=1,
shared_xaxes=True,
vertical_spacing=0.06,
subplot_titles=("Simulated returns", "Volatility clustering: squared returns vs conditional variance"),
)
fig.add_trace(go.Scatter(x=df["t"], y=df["return"], mode="lines", name="r_t"), row=1, col=1)
fig.add_trace(
go.Scatter(x=df["t"], y=df["eps_sq"], mode="lines", name=r"$\varepsilon_t^2$", line=dict(width=1)),
row=2,
col=1,
)
fig.add_trace(
go.Scatter(x=df["t"], y=df["h"], mode="lines", name=r"$h_t$", line=dict(width=2)),
row=2,
col=1,
)
fig.update_yaxes(title_text="return", row=1, col=1)
fig.update_yaxes(title_text="variance", row=2, col=1)
fig.update_xaxes(title_text="time", row=2, col=1)
fig.update_layout(height=650, legend=dict(orientation="h", yanchor="bottom", y=1.02, xanchor="left", x=0))
fig.show()
# Variance forecasts: baseline vs a one-time "shock" at the last observation
horizon = 80
eps_last = float(eps[-1])
h_last = float(h[-1])
h_fore = garch11_forecast_variance(eps_last=eps_last, h_last=h_last, omega=omega, alpha=alpha, beta=beta, horizon=horizon)
shock_multiplier = 6.0
eps_last_shocked = shock_multiplier * np.sqrt(h_last)
h_fore_shocked = garch11_forecast_variance(
eps_last=eps_last_shocked, h_last=h_last, omega=omega, alpha=alpha, beta=beta, horizon=horizon
)
h_bar = omega / (1.0 - alpha - beta)
lookback = 250
hist_x = np.arange(T - lookback, T)
fore_x = np.arange(T, T + horizon)
fig = go.Figure()
fig.add_trace(go.Scatter(x=hist_x, y=h[-lookback:], mode="lines", name=r"historical $h_t$"))
fig.add_trace(go.Scatter(x=fore_x, y=h_fore, mode="lines", name=r"forecast (baseline)"))
fig.add_trace(go.Scatter(x=fore_x, y=h_fore_shocked, mode="lines", name=r"forecast (last shock)", line=dict(dash="dot")))
fig.add_trace(
go.Scatter(
x=np.concatenate([hist_x, fore_x]),
y=np.full(hist_x.size + fore_x.size, h_bar),
mode="lines",
name=r"long-run $\bar h$",
line=dict(dash="dash"),
)
)
fig.update_layout(
title="GARCH variance forecasts: mean reversion and the impact of a market shock",
xaxis_title="time",
yaxis_title="variance",
height=500,
)
fig.show()
Notes / diagnostics#
Volatility persistence is largely controlled by \(\alpha + \beta\); values close to 1 create slow decay and strong clustering.
If \(\alpha + \beta \ge 1\), the long-run variance \(\bar h\) is not finite (often discussed as IGARCH / near-unit-root volatility).
Real data often has heavy tails; using a Student-t distribution for \(z_t\) is common in practice.